! This is a 1D advection example using square initial condition and periodic
! boundary condition for Crank-Nicolson finite difference scheme.
!
! Li Dong <dongli@lasg.iap.ac.cn>
!
! - 2018-03-20: Initial creation.
! - 2024-03-23: Use adv_1d_test_case_mod.

program crank_nicolson_adv_1d_case

  use adv_1d_test_case_mod
  use iso_c_binding

  implicit none

  interface
    type(c_ptr) function gsl_vector_alloc(n) bind(c) result(vec)
      use iso_c_binding
      integer(c_size_t), value :: n
    end function gsl_vector_alloc

    subroutine gsl_vector_set(vec, idx, val) bind(c)
      use iso_c_binding
      type(c_ptr), value :: vec
      integer(c_size_t), value :: idx
      real(c_double), value :: val
    end subroutine gsl_vector_set

    real(c_double) function gsl_vector_get(vec, idx) bind(c)
      use iso_c_binding
      type(c_ptr), value :: vec
      integer(c_size_t), value :: idx
    end function gsl_vector_get

    subroutine gsl_vector_free(vec) bind(c)
      use iso_c_binding
      type(c_ptr), value :: vec
    end subroutine gsl_vector_free

    subroutine gsl_linalg_solve_cyc_tridiag(diag, e, f, b, x) bind(c)
      use iso_c_binding
      type(c_ptr), value :: diag
      type(c_ptr), value :: e
      type(c_ptr), value :: f
      type(c_ptr), value :: b
      type(c_ptr), value :: x
    end subroutine gsl_linalg_solve_cyc_tridiag
  end interface

  real, allocatable :: rho(:,:)     ! Tracer density being advected at cell centers

  type(c_ptr) a1
  type(c_ptr) a2
  type(c_ptr) a3
  type(c_ptr) b
  type(c_ptr) y

  integer, parameter :: ns = 1      ! Stencil width
  integer i
  character(30), parameter :: scheme = 'crank_nicolson'

  namelist /params/ nx, nt, dt, u

  call get_command_argument(1, namelist_path)
  inquire(file=namelist_path, exist=is_exist)
  if (is_exist) then
    open(10, file=namelist_path)
    read(10, nml=params)
    close(10)
  end if

  allocate(rho(1-ns:nx+ns,2))
  a1 = gsl_vector_alloc(int(nx, c_size_t))
  a2 = gsl_vector_alloc(int(nx, c_size_t))
  a3 = gsl_vector_alloc(int(nx, c_size_t))
  b = gsl_vector_alloc(int(nx, c_size_t))
  y = gsl_vector_alloc(int(nx, c_size_t))

  call adv_1d_test_case_init('square', ns, rho(:,old))
  call output(scheme, 0, ns, nx, x, rho(:,old))

  ! Run integration.
  print *, time_step, sum(rho(1:nx,old))
  do while (time_step < nt)
    call crank_nicolson(rho(:,old))
    do i = 1, nx
      rho(i,new) = gsl_vector_get(y, int(i - 1, c_size_t))
    end do
    call apply_bc(ns, nx, rho(:,new))
    call advance_time()
    call output(scheme, time_step, ns, nx, x, rho(:,old))
    print *, time_step, sum(rho(1:nx,old))
  end do

  call adv_1d_test_case_final()

  deallocate(rho)
  call gsl_vector_free(a1)
  call gsl_vector_free(a2)
  call gsl_vector_free(a3)
  call gsl_vector_free(b)
  call gsl_vector_free(y)

contains

  subroutine crank_nicolson(rho)

    real, intent(in) :: rho(1-ns:nx+ns)

    real c
    integer(c_size_t) i

    c = dt / dx * 0.25d0
    do i = 1, nx
      ! Although u is constant in this case, we still write as it is variant.
      if (i == 1) then
        call gsl_vector_set(a1, i - 1, - c * u)
        call gsl_vector_set(a3, i - 1,   c * u)
      else if (i == nx) then
        call gsl_vector_set(a1, i - 1, - c * u)
        call gsl_vector_set(a3, i - 1,   c * u)
      else
        call gsl_vector_set(a1, i - 1, - c * u)
        call gsl_vector_set(a3, i - 1,   c * u)
      end if
      call gsl_vector_set(a2, i - 1, 1.0)
      call gsl_vector_set(b, i - 1, c * u * rho(i-1) + rho(i) - c * u * rho(i+1))
    end do
    call gsl_linalg_solve_cyc_tridiag(a2, a3, a1, b, y)

  end subroutine crank_nicolson

end program crank_nicolson_adv_1d_case
